
#!/usr/bin/perl
use Getopt::Std;
getopt("d",\%args);
my $dir = $args{d};
chomp($dir);
if($dir eq "")
{
        die "entered dir name where subfolders are created is empty\n";
}
open(DOS,">$dir/X_dosage");
system("gunzip $dir/23/male/*.gz");
system("gunzip $dir/23/female/*.gz");
open BUFF,"$dir/23/female/beagleout_female.geno_female.bgl.gz.gprobs" or die "can't open $dir/23/female/beagleout_female.geno_female.bgl.gz.gprobs";
open TFAM,"$dir/whole_geno_phased.tfam" or die "no file exists $dir/whole_geno_phased.tfam";
open TPED,"$dir/dosage.tped" or die " no file $dir/dosage.tped exists\n";
$i =2;
while($line = <TFAM>)
{
	chomp($line);
	@tfam = split(" ",$line);
	$tfam{$tfam[1]} = $i;
	$i = $i+3;
}	
$female = <BUFF>;
@female = split(" ",$female);
shift(@female);
shift(@female);
shift(@female);
$num = 0;
for($i=0;$i<@female;$i++)
{
	#$hash{$female[$i]} = $tfam{$female[$i]};
	$exist_dosage[$num] = $tfam{$female[$i]} -2;
	$num++;
	$exist_dosage[$num] = $tfam{$female[$i]} -1;
        $num++;
	$exist_dosage[$num] = $tfam{$female[$i]} ;
        $num++;
	$i++;
	$i++;
}
while($line = <TPED>)
{
	chomp($line);
	@tped = split(" ",$line);
	shift(@tped);
	$rsid = shift(@tped);
	shift(@tped);
	shift(@tped);
	$female = <BUFF>;
	@female = split(" ",$female);
	$rsid_fe = shift(@female);
	$allel1 = shift(@female);
	$allel2 = shift(@female);
	#print DOS $rsid_fe."\t".$allel1."\t".$allel2."\n";	
	@dosage = ();
	$dosage_num = 0;
	if($rsid_fe eq $rsid)
	{
		
		for($i=0;$i<@tped;$i++)
		{
			$A1 = $tped[$i];
			$i++;
			$A2 = $tped[$i];
			if($A1 eq "" || $A2  eq "")
			{
				die "order change while merging female x chr dosage to whole dosage\n";

			}
			if($A1 eq $allel1 && $A2 eq $allel1)
			{
				$dosage[$dosage_num] = 1;
        			$dosage_num++;
        			$dosage[$dosage_num] = 0;
        			$dosage_num++;
        			$dosage[$dosage_num] = 0;
        			$dosage_num++;

			}
			elsif($A1 eq $allel2 && $A2 eq $allel2)
                        {
                                $dosage[$dosage_num] = 0;
                                $dosage_num++;
                                $dosage[$dosage_num] = 0;
                                $dosage_num++;
                                $dosage[$dosage_num] = 1;
                                $dosage_num++;

                        }
			elsif($A1 eq "0" && $A2 eq "0")
                        {
                                $dosage[$dosage_num] = "-";
                                $dosage_num++;
                                $dosage[$dosage_num] = "-";
                                $dosage_num++;
                                $dosage[$dosage_num] = "-";
                                $dosage_num++;
				
                        }
			else
			{
				$dosage[$dosage_num] = 0;
                                $dosage_num++;
                                $dosage[$dosage_num] = 1;
                                $dosage_num++;
                                $dosage[$dosage_num] = 0;
                                $dosage_num++;

			}
		}
		for($i=0;$i<@female;$i++)
                {
			$dosage[$exist_dosage[$i]]  = $female[$i];
		}
		$dosage = join(" ",@dosage);
		print DOS "$rsid_fe $allel1 $allel2 $dosage\n";
	}
	else
	{
		die " rsid of dosage female $rsid_fe does not match with all rsid $rsid\n";
	}
}
